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Finite-time coherent sets inhibit mixing over finite times. The most expensive part of the transfer 
operator approach to detecting coherent sets is the construction of the operator itself. We present 
a numerical method based on radial basis function collocation and apply it to a recent transfer 
operator construction J8* that has been designed specifically for purely advective dynamics. The 
construction [8] is based on a “dynamic” Laplace operator and minimises the boundary size of the 
coherent sets relative to their volume. The main advantage of our new approach is a substantial 
reduction in the number of Lagrangian trajectories that need to be computed, leading to large 
speedups in the transfer operator analysis when this computation is costly. 


I. INTRODUCTION 

Finite-time coherent sets (FTCS) [9j EJ E] in non¬ 
linear, possibly time-dependent, dynamical systems are 
connected regions that are maximally dynamically dis¬ 
connected from surrounding phase space when evolved 
over a specified time duration of finite length. If the dy¬ 
namical system is an advection-diffusion equation, e.g. 
a Fokker-Planck equation, finite-time coherent sets are 
those regions for which there is minimal exchange of a 
passive tracer across their boundary. On the other hand 
if the dynamical system is a model of (purely advec¬ 
tive) fluid flow, finite-time coherent sets are those regions 
for which there is minimal exchange of fluid across their 
boundary in the presence of small noise or diffusion. The 
exchange of fluid across a boundary by small diffusion is 
proportional to the size of the boundary, and in the case 
of purely advective fluid flow, one can devise an alterna¬ 
tive characterization of FTCS which minimizes boundary 
size relative to volume [8]. 

In the pure advection setting, for volume-preserving 
dynamics, the constructions in this alternative character¬ 
ization have been shown to arise as zero-diffusion limits 
of the “classical” FTCS approach [9 . Thus, the approach 
[8] can be viewed both as an equivalent geometric reinter¬ 
pretation of the classical mathematical definitions based 
on probability and diffusion [9], and as an alternative nu¬ 
merical method for identifying FTCS. In the present pa¬ 
per we apply numerical analysis methods based on radial 
basis function collocation to very efficiently compute es¬ 
timates of the main transfer operator object [8]. Because 
of the flexibility of radial basis functions, this approach is 
suited to irregularly-shaped domains and computations 
based on trajectory data. 
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Denote by M our domain of interest at the initial 
time to. The domain M is a compact, connected sub¬ 
set of with C 1 boundary (formally, for the theo¬ 
rems quoted later, we will assume M is in fact a com¬ 
pact, connected, C°° Riemannian manifold of vanishing 
curvature). Our advective dynamics is either a nonlin¬ 
ear volume-preserving map T : M -A T{M ) (or several 
compositions of possibly different maps), or a nonlinear 
volume-preserving flow map $(x,to,t), which solves the 
ODE x = F(x,t). In the latter case, we assume that 
F(x,t) is defined on a large enough domain to advect 
our initial domain M forward for some finite time. 

In the purely advective setting, how well a set “mixes” 
is generally related to the length or irregularity of 
the boundary of the set. Measures of mixing specifi¬ 
cally designed for pure advection, such as the mix-norm 
[5j 22 , 28, penalize the geometric irregularity of advected 
passive scalar concentrations by comparing the difference 
of the concentration from a uniform concentration over 
an infinite range of spatial scales. Other advective mix¬ 
ing measures directly measure the size of the boundary 
of an initially compactly-supported passive scalar con¬ 
centration [ 15 ] . 

Our objective is to identify regions of M that have low 
boundary size relative to volume, and retain low bound¬ 
ary size relative to volume as they are advected. Such 
regions do not evolve to have highly filamented bound¬ 
aries, which under the addition of small noise would lead 
to substantial mixing; see Figure [l] By minimising the 
size of the boundary relative to volume either throughout 
the flow duration or at the beginning and the end of the 
flow, we find finite-time coherent sets, which mix least 
with the surrounding phase space. 

In recent years there have been several geometric meth¬ 
ods proposed to characterize either trajectories or co¬ 
dimension 1 surfaces that represent coherent structures 
m Ha ng nzi [2oi eh eh eu. In some two-dimensional 
cases Baozi, the aim of the geometric methods is to find 
curves that undergo small amounts of stretching. The 
new approach based on the dynamic Laplacian [8] di- 
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FIG. 1: The two-dimensional shape on the left has a 
low boundary size to volume ratio. The shapes on the 
right are three possible images of the shape on the left 
under three different nonlinear volume-preserving 
dynamical systems over a finite time duration. Under 
dynamics ‘a’ the set on the left retains a low boundary 
size to volume ratio, but under dynamics c b’ and ‘c’ the 
boundary size is significantly increased. 


rectly minimizes boundary size (in arbitrary dimensions) 
when subjected to finite time dynamics. Our main contri¬ 
bution in the present work is to reduce the computational 
effort by taking advantage of the underlying functional 
representation of the transfer operator approach. 

Our approach is similar in spirit to the work of 
Williams et aLpT! in this issue. There, the authors use 
a Galerkin scheme with globally supported radial ba¬ 
sis functions to compute the constructions in [9\. Here, 
we use collocation, locally supported functions, and the 
diffusion-free operator proposed by Froyland [8] (rather 
than the diffusive approach El) - but the goal is very 
much the same, namely constructing (data-)efficient dis¬ 
cretisations of the relevant operators in order to compute 
coherent sets. We believe that ultimately, locally sup¬ 
ported functions might yield the more efficient numerical 
scheme when aiming for high resolution. 


II. ADVECTIVE SETUP FOR FINITE-TIME 
COHERENT SETS 


ary size to volume ratios at an initial time to and a final 
time tf. 


Definition 1: Define 
h D /Ti = - ( 

{*0’ tf } ’ 2 V min {l d {M x ).l d (M 2 )} )' U 

We wish to find the minimising T: T* = 

argmin r h^ o tf |(T). The finite-time coherent set is de¬ 
fined to be whichever of Mi, M 2 has the least volume; 
the boundary is formed by the hypersurface T*. 

The above mathematical construction in the absence of 
any dynamics arises naturally in isoperimetric theory 0, 
where purely geometric questions of how to disconnect a 
compact manifold in such a way that the disconnecting 
hypersurface size relative to the volume of the smallest 
disconnected piece is minimised. Answers to such clas¬ 
sical questions reveal important information about the 
geometry of the manifold. On the other hand, in the 
present constructions the nonlinear dynamics plays the 
dominant role. We refer the reader to [8] for further fur¬ 
ther background and discussion of these ideas. 

If the flow time is long, it could be that an evolved 
region develops a filamented boundary during some in¬ 
termediate time interval, but then loses this filamentation 
by the final flow time. For the times that the filament is 
present, small-scale diffusion will decrease the coherence 
of the set. In some applications, it may be important 
to not allow such transient filaments. By extending the 
above basic characterisation of FTCS to include a series 
of times to, £ 1 ,..., tf, one can penalise such transient fil¬ 
ament at ions. Define 




iYTiZy id-i($(r,t 0 ,tj)) 

min{4(M 1 ),4(M 2 )} 


where £ n _ 1 := tf, as the natural generalisation of 

h{t 0 t f }CO* Alternatively, for continuous penalisation, we 
can consider 


h D rn _ Uo r ^-iWr,io ,t))dt 
[t °’ tf]{ ’ min{^(M 1 ),^(M 2 )} ’ 


(3) 


as a time-continuous generalisation of hj^ Q tf | (r); see |] 
for further details. 


For d > 1 we denote d— 1-dimensional volume by id-\\ 
for example, if our fluid is three-dimensional i 2 measures 
surface area, while if our fluid is two-dimensional, i\ mea¬ 
sures curve length. All of the theory and techniques out¬ 
lined here work in arbitrary finite dimensions d > 1. 

Let T C M be a smooth d — 1-dimensional object that 
separates M into two subregions Mi and M 2 . The hy¬ 
persurface T is the common boundary of Mi and M 2 , 
whose size we wish to minimise relative to the smaller 
of the volumes ^(Mi),£rf(M 2 ). Denote by <£>(•, t(M) the 
flow map from time to to time t. Below we define an 
objective function that calculates the combined bound- 


ill. A DYNAMIC LAPLACE OPERATOR 

We use eigenfunctions of a dynamic Laplace opera¬ 
tor to numerically find the minimising hypersurface T*. 
To motivate this approach, we introduce the dynamic 
Sobolev constant 18 . Define 

D m l|V/l|Li+||V(/°$(-,*f,to))||Li 
2inf aeR ||/-a|| i i 

The above constant is modelled on the Sobolev con¬ 
stant common in isoperimetric theory [3], where only 
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the first term in the numerator is present; see [8] for 
details. There is a very strong formal connection be¬ 
tween the constants infr tf }(T) (purely geometric) 

and inf/ec°°(M) s ^ 0 t f }(/) (purely functional), namely 

inf hf t0 , tf} (r) = /e mf M) sf Wf} (/) (5) 

(see [8j for the formal statement). This dynamic Federer- 
Fleming Theorem is a generalisation of the celebrated 
Federer-Fleming Theorem [7] known throughout differ¬ 
ential geometry. 

Our numerical approach will be to find an / minimis¬ 
ing tf j(/). In general, the infimum in the RHS of (|5| 
is only approached by a sequence of C°° functions, and 
there is no simple formula for such a sequence of func¬ 
tions. On the other hand if we use the L 2 norm instead 
of the L 1 norm in the RHS of ([5|, the infimum is attained 
by a smooth function, which is the second eigenfunction 
of the dynamic Laplace operator, defined by 

A d := (A + P*AP)/2. (6) 

In ([6]), A is the standard Laplace operator, Vf = 
/ o <£>(•, £f, to) is the Perron-Frobenius (or transfer) op¬ 
erator for <£>(•, to,£f), and V*f = / o <!>(•, to,£f) is its 
dual, the Koopman operator for <!>(•,to,tf). Note that 
the first Laplace operator in © acts on M, while the 
second Laplace operator acts on <I>(M, to,tf); this is 
crucial for A D to feel the geometry at both time in¬ 
stants. In fact, V*AV is the natural pullback of the 
Laplace operator on <F(M,to,tf) under <!>(•,to,tf), and 
is the Laplace-Beltrami operator on M with respect to 
the pullback of the Euclidean metric on <F(M, to, tf). Let 
S denote the trivial (Euclidean) Riemannian metric on 
<F(M, to, tf); pulling S back under <£>(•, to, tf) we obtain 
the Riemannian metric $(•, to, tf)*5 on M, and the map 
$(Vo,if) : -> (,®{M,t 0 ,t f )),5) is an 

isometry. Then A 4>( . )toitf) . 5 / = (A s (f o $(•, t 0 , t f )) _1 ) o 
<£>(•, to, tf) = V*A$T , where As is the Laplace-Beltrami 
operator on the Riemannian manifold (<F(M,to,tf)),5) 
see Section 3.2 [8] or p27 [2]. 

If M has a boundary, then the eigenfunction f : M 
M is required to satisfy generalised (oblique) Neumann 
boundary conditions. Denote by n(x) the unit outward 
normal at x G dM. We require 

V/(*) • [(/ + )n(aO] =0, x e dM, (7) 

where C Xtto ,t r ■= D$(x,t 0 ,tf) T ■ D®(x,t 0 ,t f ) is the 
(right) Cauchy-Green deformation tensor for the trans¬ 
formation <£>(•,to,tf). This boundary condition also has 
a natural pullback interpretation, see Section 3.2 [9,. 

A. Spectral properties of the dynamical Laplacian 

The spectral properties of the L 2 eigenproblem ©-0 
are (see f8j for formal statements): 


1. A D is self-adjoint. 

2. The eigenvalues form a decreasing sequence 0 = 
Ai > A 2 > ■ ■ • with A n —)► —oc. 

3. The leading eigenfunction fi = 1, and the eigen¬ 
functions /i, / 2 ,... corresponding to distinct eigen¬ 
values are pairwise orthogonal in L 2 . 

4. One has the following variational characterisation 
of eigenvalues: if /i, / 2 ,... are arranged to be or¬ 
thonormal, denoting Xk = span{/i, / 2 , 

l|V/||| 2 + ||V(/o$(.,t f ,t 0 ))||2 2 

m\\h 

:{f,fi)=0,i = l,...,k-l.} (8) 

with the infimum achieved only when / = /&. Note 
that setting k m 2 minimises © when || • || L i is 
replaced with || • \\ L 2 . 

From (|8|, one sees eigenvalues that are far from 
zero correspond to functions / that have high gradient 
over a large part of the domain (we call such functions 
“irregular” as they are far from the most regular function, 
namely the constant function). The equality © makes 
this exact for k = 2, and in fact, one also has a dynamic 
Cheeger inequality 0: 

in f hf to , f} ( r )<2v^, (9) 

which links the least possible dynamic boundary size to 
volume ratio with the first nontrivial eigenvalue A 2 of the 
dynamic Laplace operator A D . 

If there are K > 1 independent ways to disconnect M, 
such that the disconnections each have similarly small 
dynamic boundary size to volume ratio, then there should 
be a cluster of eigenvalues A 2 ,..., \k+i that are much 
closer to 0 than A^+ 2 > A^+ 3 , • • •• In such a situation, the 
K finite-time coherent sets can be extracted from level 
sets of the eigenfunctions /&, fc = 2,...,AT + l. Perhaps 
the simplest way to do this is to find the K level sets T 
by optimising h D (r) with K separate line searches. For 
related approaches in the autonomous dynamics setting, 
see also mm- 

B. Coherent sets and objectivity 

Existing transfer operator methods for identifying 
finite-time coherent sets o nu na ns extract the co¬ 
herent sets from super- or sub-level sets of singular vec¬ 
tors of the transfer operator corresponding to singular 
values close to 1. We follow the same strategy here, us¬ 
ing level sets of the solutions of the eigenproblem ©-© 
that correspond to eigenvalues close to 0. The numerical 
algorithm is detailed in Section 

Concerning objectivity of the algorithm, it is shown 
in [8] that if the domain M is additionally subjected to 
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a time-dependent affine transformation with orthogonal 
linear part, the solutions of the corresponding eigenprob- 
lem ©-© are simply transformed versions of the solu¬ 
tions of the original eigenproblem. Thus, as the coherent 
sets are extracted from level sets of the eigenfunctions of 
©- 0 , the extracted coherent sets are also identically 
transformed as required for objectivity of the method. 


C. Penalising boundary size at multiple times 


If one wishes to penalise the boundary size at several 
discrete times between £ 0 and t f then define 


n— 1 




( 10 ) 


i=0 


where Pt 0 ,nf = /°$(-, U, to)- Consider the eigenproblem 
A£, tl ,... ltn _ a ,t f /(*) = A f(x), x € M, (11) 
with boundary condition 


V/(x) • 


"n—1 




.i=0 


= 0, x £ dM. 


( 12 ) 


If one wishes to penalise the boundary size at all times 
between to and t f then define 

f Q K,AV t0 , t dt, (13) 

where 'Pto, tf = f o 4>(-, C to). Consider the eigenproblem 

A \to,t t ]f( x ) = A f(x), xeM, (14) 

with boundary condition 


V/(x) • f C x l t n(x) 
IJo 


dt 


= 0, x G dM. (15) 


In both of the above cases, the obvious versions of © 
and © hold [8]. 


IV. DISCRETIZATION OF THE DYNAMIC 
LAPLACIAN EIGENPROBLEM 


Following Platte and Driscoll [26 , we are going to ap¬ 
proximate the eigenproblem ©H [7]) by collocation with 
radial basis functions. Here, we choose the Wendland 
functions ^ = 2pd,k • [0, 00 ) M for various k G N, which 
have support on [0,1], are polynomials of a certain de¬ 
gree and lead to strictly positive definite interpolation 
matrices [29]. 


A. The case 4>(M, to,tf) = M 

We first choose a set Y = fei,..., 2 /n}cMof centers. 
The corresponding basis functions <pj : M M are given 
by 


<Pj(x) = i’(4 x ~yjh), 

j = 1,..., n, where e > 0 is the shape parameter which 
scales the size of the support. 

The eigenfunctions / of A D will be approximated by 
functions from A = span(<^i,..., ip n ) C L 2 (M), i.e. by 
linear combinations of the form 

n 

K x ) = '52<xj<Pj( x )- 

i=1 

To this end, we choose Z inner collocation nodes X{ n = 
{xi,..., aq} C M in the interior of M as well as m = n—£ 
nodes Xbd = { 2 ^+ 1 ,... ,x n } on the boundary of M. If 
M has no boundary, then £ = n and Xbd = 0- Now, 
given a (sufficiently smooth) function / : M —>■ M which 
satisfies the boundary condition 

£bd/ := V/(:e) • [(/ + )n(x)] = 0 

on 3M, the coefficients = (oq,..., <a n ) T of its interpo¬ 
lating approximation f £ A are given by the solution of 
the linear system 

A(X = Eofin 

where 



Ai„ 


I 


~f( x l)" 

A = 

Fbd 

j Eq = 

0 

, /in = 



and 

^in = {<Pj i = 1,...,£, j = l,...,n 

Tbd = (Z^bd^jf (^i))^j 5 1,...,77/, J = l,. 

For the case 9M = 0, the fact that the interpolation 
matrix Ai n is (strictly) positive definite implies the in¬ 
vert ibility of A. In the case dM ^ 0, however, invert- 
ibility does not seem to be guaranteed for an arbitrary 
choice of the nodes Aq n ,Xbd [26] . In our numerical ex¬ 
periments, however, A was always invertible. 

Given some arbitrary linear operator C : L 2 (M) 
L 2 (M), the image Cf of a function / G A is given by 

n 

£f( x ) = '52<xj£<Pj( x )- 

3 = 1 

Correspondingly, the values gi = Cf(xi) of this image 
at the inner collocation nodes aq, i = 1 ,...,£, can be 
computed by 

9 = -Zin^ = Li n A Eofin, 
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where 

L- m = (Cipjixi))^ , i = 1,... j = 1,... ,n. 

That is, the matrix L := Li n A~ 1 E 0 can be seen as a 
discrete version of the given linear operator C. 

In our case, we are concerned with the linear operators 
V, A and V * (cf. [6| and so the approximating operators 
are P = P[ n A~ 1 Eo with 

^Pin = ('P<Pj(xi)) ij , i = 1,... j = 1,.. • ,n, 

P* = P^* n A _1 P 0 with 

and D = Di n A _1 Po with 

Pin = i = 1,..., £, j = 

The discrete eigenproblem associated to ©-0 is 

l( D + P* DP )f = Xf. 

We note that P is unitary, so P* = P _1 . In general, 
this property is not inherited by the matrix P which is 
constructed via collocating the action of T _1 . Likewise, 
the matrix P* is constructed by collocating the action of 
T, i.e. one may view P* as taking the dual of P followed 
by collocation. In general, this approach will not guar¬ 
antee that P* is the dual of P. As a result, the discrete 
dynamic Laplacian ^(D + P*DP) is not necessarily self- 
adjoint and therefore its spectrum not necessarily real (in 
contrast to A d )[2S]. 

However, if we interchange these two operations, i.e. if 
we first collocate and then take the dual, we obtain P T 
which is the dual of P. In our numerical experiments, us¬ 
ing P T instead of P* turned out to yield a much improved 
spectrum, i.e. in Section [V] we solve the eigenproblem 

1 -{D + P T DP)f = \f. (16) 

Remark. Note that, as mentioned above, this dis¬ 
cretization applies to an arbitrary linear operator C (cf. 
Platte and Driscoll [26]). In particular, one may apply 
this to the “classical” approach for computing coherent 
pairs [9], namely discretizing the transfer operator of a 
small random perturbation of the underlying system. In 
this case, in addition to the matrix P, one would com¬ 
pute a discretization D £ of some local averaging operator 
V s and compute the largest singular values and vectors 
of the matrix D £ PD £ . 


made which strongly influence the resulting approxima¬ 
tion quality. Here, we collect a few comments on a con¬ 
crete implementation. Corresponding Matlab codes can 
be obtained from the authors. 

We first need to decide on a set Y of centers for the 
basis functions. While in general, scattered points are 
fine for approximations with radial basis functions (in 
fact, this is often one of the main motivations for using 
radial basis functions), in our experiments regular grids 
led to a much improved accuracy. In our 2D examples, 
a set Y of 1000-2000 centers led to sufficiently accurate 
results. 

The next choice concerns the shape parameter 5 which 
governs the radius of the support of the basis functions. 
Again, the computational results are typically highly sen¬ 
sitive with respect to this parameter. In general, better 
results are obtained when choosing a smaller value for 
e , leading to larger supports. On the other hand, the 
condition number of the interpolation matrices (i.e. Ai n , 
etc.) increases with decreasing 6 and for values too small, 
instabilities will occur. In our experiments, we chose e 
such that the support of a basis function overlapped with 
roughly 100-200 other ones. 

For the collocation points Xi n , Xbd, the same remarks 
apply as for the centers. Here, we use choose Xi n UAbd = 
Y since this seems to preserve the spectral properties of 
the various matrices best. Only in the second example 
(cf. Section [V]), we slightly shifted the centers away from 
the boundary in order to keep the code simpler. Note 
that the choice X{ n U Xbd = Y requires a special treat¬ 
ment of those entries of the discrete Laplacian where the 
collocation point coincides with the center (because of 
the zero in the denominator of the Laplacian in polar co¬ 
ordinates, which we employ due to the radial structure 
of our basis functions). 

To compute the inverse of the Cauchy-Green ten¬ 
sor, we integrate the variational equation backward in 
time, i.e. explicitly compute D x $($(x,toAf)AfAo) f° r 
x G W>d backwards along the nonlinear trajectory 
{<F(x,£o P)}t 0 <t<t f computed prior in forward time. For 
longer flow times, this may yield a (nearly) singular ma¬ 
trix. 

We also have to compute the various matrices: For 
Abd and Pi n , we first apply the differential operators to 
the basis functions x c pj(pc ) and then evaluate at the 
collocation points. Similarly, when computing Pi n , we 
compute 

(-Pin)zj = f,to))- 


B. Implement at ional details 

In implementing the approach described above, var¬ 
ious choices on the numerical parameters have to be 


Finally, for the resulting discrete eigenproblem (16), 
we compute a few eigenvalues of smallest magnitude by 
using an implicitly restarted Arnoldi method [19] as im¬ 
plemented in the eigs function in Matlab (i.e. with pa¬ 
rameter ; SM ; ). 
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C. The general case <h(M, to,£f) p M 

Whenever M := <L(M, £ 0 pf) <f_ M, it will be necessary 
to work with two distinct approximation spaces that are 
supported on M and M, respectively. In the context of 
radial basis function approximations, it appears natural 
to define the set of centers Y in M by the images of the 
centers in Y under the flow (in this section for brevity 
we will write for <!>(•, to, tf)) which here we assume to be 
a bijection (as is the case if the underlying vector field F 
is locally Lipschitz w.r.t. x). In fact, the discrete dynamic 
Laplacian attains a particular simple form in this case as 
we will see now. Choosing Y = {&(y) \ y E V}, we have 
the associated basis functions 

v>j( x ) = ^(4 x -yjh), 

j = 1 ,...,n, so that we can approximate functions 
/ : M —>■ M by elements from the space A = 

span(<£i,..., (p n ). Similarly, we define the set of col¬ 
location points as X- in = <h(Xi n ) = £ n } = 

{^>(xi),..., $(x n )} and X hd = <$>(X hd ). 

For the case dM = 0 one can proceed as follows: Given 
some function f{pc) = Y^j=i a j ( Pj( x ) £ A, the approxi¬ 
mation of Vf within A will be given by a linear combi¬ 
nation of the form YAj=i I n fact, requiring this 

equality to hold in the collocation points £*, we obtain 
the system 


n n 

3 = 1 3 =1 

i = 1 ,...,n. Since ~ 1 (xi ) = ^> _1 (^>(^)) = Xi, we 

obtain the matrix representation 

p a = i- 1 ^, 


with D' = (A(pj(x)\ x= ^ i )ij and the discrete dynamic 
Laplacian is 


V D ° 


P*D a P a ). 


(17) 


Note that the matrix 0 is a mapping on the coef¬ 
ficient vectors a of some RBF-approximation / E A - 
in contrast to (16), where the matrices are mappings on 
vectors f[ n G of function values. Certainly, we can 
rephrase the eigenproblem here in terms of /i n , too. To 
this end, we apply A and A~ x (resp. A and A -1 ) appro¬ 
priately, yielding the matrices 


P = AA^AA- 1 = I 
p* = aa- 1 !!- 1 = I 
D = AA^D'A- 1 = D'A- 1 


D = AA^d'A - 1 = D'A~\ 


i.e. we obtain the discrete dynamic Laplacian 

\{D + D). (18) 


Let us explore a little the similarity between 
pressions § and Jl8| ). In the second term of 


the ex- 
(18l) we 


have (D') i:j = A(pj(x)\ x=Ai = A^ j {\\x-^{y j )h)\x=^{x i )^ 
That is, in computing P', we are using the distance ma¬ 
trix (||^>(x^) — &(yj)\\ 2 )ij of pairwise distances between 
the images of the centers and the collocation points. This 
is consistent with the fact that the operator P* AV is the 
Laplace operator on M endowed with the ^-pullback of 
the Euclidean metric on <L(M); see Section 3.2 [8j or p27 
0 . 


D. Extraction of coherent sets from trajectories 


for the discretization of V as a map (on the coefficients a) 
from A to A, where Aij = (< pj[xi )) and Aij = (<£j(£*)), 
i, j = 1,..., n. Similarly, for P*, we get the matrix 

P* = A _1 A. 

Again, one might want to use Pj instead of P* in or¬ 
der to enforce self adjointness of the discrete dynamic 
Laplacian, cf. the comment at the end of Section |IV A| 
For the application of the Laplacian A, we again have 
to distinguish whether it is being applied to a function 
from A or A. Mimicking the above derivation for P, we 
obtain the matrix representation 

D a = A~ 1 D' 

for the discretization of A as a map (again, on the co¬ 
efficients a) from A to A , where D' = (A (Pj(x)\ x=x .)ij. 
Likewise, for A (we call the second A in ^ A as it acts 
on $(M)) 


D a = A- 1 !)' 


The algorithm we use in the following two-dimensional 
volume-preserving case studies is described below. We 
consider the case <I>(M, £oPf) = M and correspondingly 
use the approach described in Section |IV A| Algorithm 
1 is presented for the situation where the boundary size 
is compared at the initial and final times; the obvious 
extensions can be made if additional comparison times 
in the interval [to At] are desired. We outline the main 
steps here, more detailed comments on the individual 


steps can be found in Section IV B 


Algorithm 1: 

1. Choose a set of centres Y = {yj}™^ in M and 
shape parameter e > 0. 


2. Select a set of internal collocation points X 1Y1 = 

{xi}\ =1 and a set of boundary collocation points 
.Abd = {xi}i =£+1 i n If M has no boundary, 

then £ = n and Xbd = 0. 

3. If Xbd A 0? calculate the unit normal n (xA and the 
Cauchy-Green tensor C Xi) t 0) t f for i = £ + 1,..., n. 
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4. Form the matrix P for the transfer operator and the 
matri x D for the Laplacian as described in Section 
as well as D' = (D + P T DP)/ 2. 


IV A 


5. Solve the eigenproblem ( p~ 6 ] ) ? i.e. compute several 
eigenvalues Ai > A 2 > • • • of smallest magnitude 
(e.g. using Matlab’s eigs command with the option 
’ SM ; ) and associated eigenvectors / 1 , / 2 , • • • of D'. 


6 . Iteratively scan over values of from min /2 to 
max/ 2 . For each value, extract a level curve T in 
M using Matlab’s contour function (this function 
returns a collection of points representing corners 
of a polygonal curve). To compute 4>(r, £ 0 , £f), 


(a) Either: map the points representing T directly 
with <F(-,£ 0 ,£f), 

(b) Or: compute P /2 and extract <f>(r,£o,£f) us¬ 
ing Matlab’s contour function with the same 
level set value as for T. 


7. Optimise h D (r) by running over all curves T 
formed from level sets of in Step 6 . The length 
of T and <F(r,£o,£f) are computed as the lengths 
of the polygonal curves comprising them. Report 
the T and <F(r,£o,£f) that yield the lowest value of 

h D (r). 


V. NUMERICAL EXPERIMENTS 


A. Standard map on the torus 


We consider the standard map 

T(x, y) = (x + y + a sinx, y + a sinx) (mod 27r), (19) 

with parameter a = 0.971635, a value above which a 
prominent KAM curve is destroyed, and the phase space 
shows both regular and chaotic motion; see Fig. [ 2 ] We 


27T 


^5 


2n 


FIG. 2 : Orbits of the map (19) showing regular and 
chaotic regions. 



namics of two iterates of T. The geometric action of T 2 
in phase space is indicated in Figure [3] 



0 27r 0 27r 

x x 

FIG. 3: Left: An initial colouring of phase space. Right: 
The image of the colouring under T 2 . 


We applied Algorithm 1 as follows: 1. We choose the 
centers Y = hZ 2 fl[0, 27r] 2 on a regular grid with h = 0.33, 
leading to 400 points in V, use the Wendland function 
ip 6,4 f° r the basis functions and £ = 0.4 for the value of 
the shape parameter, such that the support of each basis 
function overlaps with roughly 200 other basis functions. 

2 . We choose X[ n = Y for the inner collocation points 
and since dM = 0 we have Xbd = 0 and we do not 
need to compute any normals or Cauchy-Green tensors. 

3. We compute the matrices P and D as described in 
Sections |IV A| and |IVB| and 4. solve the eigenproblem 

Jl6). 


The resulting four leading eigenvalues were (to three 
significant figures): Ai = 6 • 10 -5 , A 2 = —1.15, A 3 = 
— 1.17 and A 4 = —2.10 which appear to be correct when 
compared to the results of a highly accurate computa¬ 
tion using a spectral method. Fig. 0] shows the spectrum 
of the discrete dynamic Laplacian 1^(D + P T DP). Note 
that this spectrum is not real (as it should be). However, 
since the part of the eigenspectrum we are interested in is 
real, namely in a neighborhood of the origin (see Fig. [4] 
(right)), we do not need to employ further techniques 
to reduce the imaginary components of the spectrum of 
^(D+P T DP) (such as using a least squares approach in¬ 
stead of interpolation|271). One could also enforce a real 
spectrum by using (D + D T )/2 instead of D. The ques¬ 
tion of how these two techniques change the spectrum is 
currently under investigation. 


•1(T 3 
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are interested in finite-time dynamics rather than time- 
asymptotic dynamics and we choose to analyse the dy- 


FIG. 4: Standard map: Spectrum of \(D + P T DP) 
(left: entire spectrum, right: eigenvalues closest to 0 ) 




























Note that in order to obtain this accuracy, we only 
evaluated the map 20-20 = 400 times. Experimentally, in 
order to obtain the same accuracy with Ulam’s method 
one needs to use a grid of 64 x 64 = 4096 boxes and 
7 x 7 = 49 sampling points per box at least, requiring 
the mapping of 4096 • 49 ~ 2 • 10 5 points. On the other 
hand for our current approach, setting up the distance 
matrices as described in Section |IV A| requires additional 
computational effort and the resulting eigenproblem is 
not sparse, in contrast to the one resulting from Ulam’s 
method. 

The eigenvalues A 2 and A 3 are of similar value, indi¬ 
cating that there are two independent ways to discon¬ 
nect M such that the disconnections each have similarly 
small dynamic boundary size; see Figs. [5] and [6j The 



0 2tt 0 27 r 


FIG. 5: Two iterations of the standard map: 
Eigenvectors (left) and fs (right) of (16), using 
n = 400 centers on a regular grid with 6 = 0.4. 


first two nontrivial eigenvectors are shown in Figure |5j 
In producing this figure (as well as all subsequent figures 
of eigenvectors of the dynamic Laplacian), we evaluated 
the RBF representation of the eigenvectors on a grid of 
100 x 100 points. Their images P /2 and P /3 under P 
are shown in Figure [ 6 j these vectors correspond to the 
geometry of the phase space after two iterations of T. 



0 2ix 0 2tt 


FIG. 6: Images P/2 (left) and P/3 (right) of the 
eigenvectors from Figure [ 5 ] 


Step 6 of Algorithm 1 involves scanning over the level 
sets of /2 (shown in the left in Figures [ 5 ] and [6| and 
calculating the length of the level set and the length of 
its image. We take the approach of version (b) in Step 
6, and Figure [7] shows the value of 0 as a function of 
7 £ M, where T 7 = {x G T 2 : fz(x) = 7 }. 



FIG. 7: Standard map: Graph of the dynamic Cheeger 
constant ([!]) in dependence of the level set value 7 . 


The minimising 7 = —0.0115 was chosen by scanning 
through 1000 equally separated values of 7 from min /2 
to max/ 2 , where is scaled so that \f 2\00 = 1- The 
curves along the corresponding level sets in and P /2 
are shown in Figure [ 8 j having lengths G(T 7 ) = 14.90 
and ^i(4>(T 7 ,to,tf)) = 13.54, while ^(Mi) = 20.12 for 
the corresponding finite-time coherent sets M\ shown in 
Figure^! so that h^ 2 (T 7 ) = 0.70. Note that the curves in 
Figure |§] (left) and (right) are similar, but not identical. 
There is no prescription for the curves to be invariant 
under T 2 ; the fact that they are approximately invariant 
is likely tied to the asymptotic dynamics of T. 




0 2tt 


FIG. 8 : Minimising curves before (left) and after (right) 
the application of T 2 . 
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The curves in Figure [8] are in the vicinity of the “in¬ 
ner” boundaries (the boundaries closer to y = 0.5) of the 
chaotic regions in Figure [2] The red features in Figure [8] 
also correspond reasonably closely to the largest regular 
regions in the upper and lower parts of Figure [2] Figure [2] 
displays the time-asymptotic dynamics of T, while we are 
analysing two iterations of T. If we increased the num¬ 
ber of iterations in our analysis, it is likely that the level 
sets of / 2 will more closely reflect the time-asymptotic 
dynamics. 

The third eigenfunction / 3 shows further finite-time 
structures that undergo low levels of deformation (the 
boundaries of the approximately red and blue sets remain 
small at the initial (Figure [5] (right)) and final (Figure [6] 
(right)) times, respectively. The total boundary lengths 
of the red/blue interfaces in Figures [5] and [6] are larger 
for / 3 than for / 2 (albeit only slightly so). These finite¬ 
time structures are highly mobile from time t = 0 to time 
t = 2 and are not at all evident in Figure [2| 

Order of convergence. One of the benefits of using 
radial basis functions for function approximation is that, 
depending on the “basic” function if one chooses, high 
convergence orders can be obtained. Here, we have cho¬ 
sen the compactly supported Wendland functions if = 


rors. With if = -06,4 we obtain an error of approximately 
10 -3 for 1 /s & 1.6. Recall that we only needed to evalu¬ 
ate the map n = 400 times here. On the other hand, for 
if = ^ 6,4 and if = ^5,3, a further decrease in 5 does not 
yield a smaller error. 



FIG. 11: Standard map: Dependence of the maximal 
absolute error in the first 4 eigenvalues on the mesh 
widths for e = 0.8. Green: if 3 , 1 , black: if 4 , 2 , red: ^ 5 , 3 ? 
blue: if 6 , 4 ). 



FIG. 10: Standard map: Dependence of the maximal 
absolute error in the eigenvalues Ai,..., A 4 on the 
support radius of the basis functions for n = 400 centers 
on an equidistant grid. Green: if 3 , 1 , black: if 4 , 2 , red: 
ip 5 , 3 , blue: ip 6A ). 


We next analyse how the eigenvalue error depends on 
the number of centers, i.e. the mesh width. To this end, 
we need to fix the value of the shape parameter 5 since it 
is known that an interpolating approximation with RBFs 
in which one scales the support of the basis functions pro¬ 
portionally to the mesh width does not converge|6j f30] . 
We perform an experiment with s = 0.8 and show the 
maximal error in the four eigenvalues with smallest mag¬ 
nitude in Figure mi We observe convergence orders of 
0 (s~ 3 ), 0 (s~ 3 ' 5 ), 0 (s~ 6 ' 5 ) and (D(s~ 9 ) as indicated by 
the black lines (from top to bottom) by using ^ 3 , 1 , ^ 4 , 2 , 
^ 5,3 and if 6 , 4 , respectively. Clearly, using basis functions 
of higher smoothness pays off, but surprisingly, if 3,1 and 
if 4,2 seem to converge with approximately the same rate 
(albeit ^ 4,2 delivers a smaller error). Finally, in Fig- 
the error for ^ 6,4 decays down to 8 • 10 -4 with 
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decreasing mesh width, but then starts to rise again for 
some unknown reason. 


B. Cylinder flow 


We first investigate how the absolute errors in the 
largest few eigenvalues depend on the shape parameter. 
Figure [lO] shows the results of an experiment for n = 400 
centers. Here, the results of the RBF approach are com¬ 
pared to those of a computation using spectral colloca¬ 
tion which appears to be highly accurate. In this figure, 
the maximal error in the four eigenvalues with smallest 
magnitude in dependence of the radius 1 /s of the support 
of the basis functions <fj is shown for various choices of 
the basic function ip Sj k- Clearly, these errors sensitively 
depend on e and smoother basis functions yield lower er- 


As a second example, we reconsider a genuinely nonau- 
tonomous system with nonempty boundary [ 12 ]: a flow on 
the cylinder M = S 1 x [0, 7 r], given by the vector field 

x = c — Aft) sin(x — vt) cos fy) + sG(g(x , y : t)) sin(t/ 2 ) 
y = Aft) cos(x — vt) sin fy) 

with Aft) = 1 + 0.125sin(2\/5t), Gfif) = 1 /{if 2 + l) 2 , 
g(x, y , t) = sin(x — vt) sin fy) + y/2 — 7r/4 and parameter 
values c = 0.5, v — 0.25 and 5 = 0.25. We consider the 
flow map T := <!>(•, to, tf) with to = 0 and tf = 40 which 
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we approximate by 400 steps of the classical Runge-Kutta 
scheme of 4th order. 

The flow exhibits strong mixing apart from two eddy¬ 
like structures which roughly retain a constant y value 
while moving along the circle (i.e. the x~) direction. Fig- 
(left) shows one of these two coherent sets as 


ure 
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a black curve, together with its image under the flow 
(right). On the other hand, selecting a set 
away from the objects identified in Fig. [l7j such as in 
Fig. I2l reveals chaotic motion. There, the set within the 


black curve from Fig. El is shifted by 2 in x-direction at 
to = 0, and shown (Fig. El left) together with its image 
at tf = 40 (right). Evidently, this pair is not coherent. 



FIG. 12 : Cylinder flow: Non-coherent pair. The set to 
the left is the one from Fig. [l7] (left), shifted by +2 in 
x-direction. 


We again see that A 2 and A 3 are relatively close, com¬ 
pared to lower eigenvalues. This indicates the presence of 
two distinct finite-time coherent sets. The corresponding 
eigenfunctions and fs are shown in Figure [l4j while 
Figure [15] shows their push-forwards P /2 and P/ 3 . 

7r 

53S 


0 ™ 27T 0 ™ 271- 




FIG. 14: Cylinder flow map: Eigenvectors /2 (left) and 
fs (right) of (16), using n = 2500 centers on a regular 
grid with e = 2 . 



We applied Algorithm 1 with the following parameters: 
1. For S = 10 -6 , we choose Y = Y 6 = ((h x Z+5) x h y Z)fl 
(.S 1 x [0, 7 r]) with h s x = (27r — 2<5)/50 and h y = 7t/ 50 
for the centers, i.e. 2500 points on a regular grid which 
is slightly shifted away from the boundary, and s = 2 
for the shape parameter such that the supports of the 
basis functions overlap with roughly 100 others. 2 . We 
choose X in = Y°\{(x,y) \ y = 0 or y = tt} for the inner 
collocation points and Xbd = T°\{(x, y) \ y 7 ^ 0 and y 7 ^ 
7 r} for the boundary points. The Cauchy-Green tensor 
was computed by integrating the (backward) variational 
equation for each point in Xbd- 3. We compute the 
matrices P and D as described in Sections llV Al and IIVBI 
and 4. solve the eigenproblem ( fl 6 | ). 

The resulting four eigenvalues with smallest magni¬ 
tude were Ai = 0.25, A 2 = —5.90, A 3 = —8.10 and 
A 4 = —22.4. Fig. p~3| shows the spectrum of the discrete 
dynamic Laplacian^ (P + P T DP). Again, the spectrum 
is not real, but it is real close to 0 which is the part we 
are interested in (cf. our comments on the spectrum in 
Section |V Ah . 
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FIG. 13: Cylinder flow: Spectrum of |(P + P T PP) 
(left: entire spectrum, right: eigenvalues closest to 0 ). 
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Figure 14 


Again, as described in Step 6 of Algorithm 1 , we scan 
over the level sets of (shown to the left in Figure [l4|) 
and compute the length of the level set and the length of 
its image. To illustrate the first of the two options in Step 
6 of Algorithm 1 , we now use Step 6 (a) and directly map 
points found on the curves T 7 = {x E S 1 x [ 0 , 7 r] : h (x) = 
7 }, 7 G [—1,1], to evaluate hj^ 40 |(r 7 ). Figure fl 6 | shows 
the value of 0 as a function of 7 . The minimising 



-1 -0.5 0 0.5 1 

7 


FIG. 16: Cylinder flow: Graph of the dynamic Cheeger 
constant (|T|) in dependence of the level set value 7 . 

7 = 0.4372 was chosen by scanning through 100 equally 
separated values of 7 from min /2 to max/ 2 , where /2 
is scaled so that \f 2\00 = 1- The curve along the cor¬ 
responding level set in /2 and its image under <!>(•, 0,40) 
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are shown in Figure [TF| having lengths G(r 7 ) = 4.52 and 
£i(<h(r 7 ,£ 0 ,tf)) = 5.01, while ^(Afj.) = 1-57 for the cor¬ 
responding finite-time coherent set M\ (red in the Figure) 
so that h^ 40 (r 7 ) = 3.04. 




FIG. 17: Minimising curves before (left) and after (right) 
the application of <£>(•, £q, tf*). 


One could now threshold fo according to Steps 6 and 
7 of Algorithm 1 to find an optimal T from f%. This 
would yield the union of the curve already found from f 2 
and a closed curve surrounding the right-hand blue set 
in Figure [T7| 


VI. CONCLUSION AND FUTURE 
DIRECTIONS 

We presented a fast method for the discretization of 
transfer operators based on collocation with radial basis 
functions (RBFs), and applied this method to compute 
finite-time coherent sets using the new advection-only 


construction from [8]. In particular, by choosing suffi¬ 
ciently smooth kernels, we observed that very high con¬ 
vergence orders can be achieved. These rapid construc¬ 
tions of accurate numerical approximations of transfer 
operators alleviate the major computational expense in 
algorithms to detect finite-time coherent sets. 

We demonstrated that we could construct accurate es¬ 
timates of the boundaries of finite-time coherent sets in a 
nonlinear cylinder flow over a rather long time duration 
using only a 50 x 50 grid of initial and final points, and 
without exploiting any dynamics or geometry of the flow 
in our choice of RBF centres or radii. Beyond construct¬ 
ing a discretisation of the main boundary value problem, 
our numerical contributions include (i) stable methods to 
handle the boundary conditions that avoid matrix inver¬ 
sion, and (ii) methods to deal with the situation where 
the entire domain evolves under the dynamics. 

In our numerical experiments we observed that the ac¬ 
curacy of the RBF collocation procedure can be sensitive 
to the various parameters involved, i.e. to the choice of 
the centers, the collocation points, and the shape param¬ 
eter 6, which controls the radius of the supports of the 
basis functions. In practice, it is therefore advisable to 
repeat a particular computation with different values for 
the parameters. We observed that choosing the colloca¬ 
tion points different from the centers does not yield a real 
spectrum of the discrete Laplacian. A more systematic 
study of this phenomenon will be the subject of future 
work. 
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